


clear
capture log close
set more off

global path ""
cd  ""
global dta1 ""

******* cpi 2
. clear

. use "bls_prices_2.dta"

* normalize cpi to equal one in base years

foreach num in 1990  {
	gen norm`num' = cpi if year==`num'
	sort year
	replace norm`num' = norm`num'[_n-1] if norm`num'==.
	gsort -year
	replace norm`num' = norm`num'[_n-1] if norm`num'==.
	gen cpi_`num' = cpi/norm`num'
	
}
* save
keep year cpi_*
sort year
save cpi_allyears_2, replace 

clear
use cps_00003
sort year serial pernum
merge m:1 year using cpi_allyears_2
*keep if _merge==3
drop _merge

* calculate real income

gen hhincome_1990 = hhincome/cpi_1990

** keep household head only
keep if relate==101
	
* Census division
gen division = 1 if region==11
replace division = 2 if region==12
replace division = 3 if region==21
replace division = 4 if region==22
replace division = 5 if region==31
replace division = 6 if region==32
replace division = 7 if region==33
replace division = 8 if region==41
replace division = 9 if region==42

save income_b1990.dta,replace


clear
use income_b1990.dta

forvalues i=1990/2020 {
use income_b1990.dta

* create real income groups

	gen hhfaminc_1990 = 1 if hhincome_1990<3000
	replace hhfaminc_1990 = 2 if hhincome_1990>=3000 & hhincome_1990<4000
	replace hhfaminc_1990 = 3 if hhincome_1990>=4000 & hhincome_1990<5000
	replace hhfaminc_1990 = 4 if hhincome_1990>=5000 & hhincome_1990<6000
	replace hhfaminc_1990 = 5 if hhincome_1990>=6000 & hhincome_1990<7500
	replace hhfaminc_1990 = 7 if hhincome_1990>=7500 & hhincome_1990<9000
	replace hhfaminc_1990 = 8 if hhincome_1990>=9000 & hhincome_1990<10000
	replace hhfaminc_1990 = 9 if hhincome_1990>=10000 & hhincome_1990<11000
	replace hhfaminc_1990 = 10 if hhincome_1990>=11000 & hhincome_1990<12500
	replace hhfaminc_1990 = 12 if hhincome_1990>=12500 & hhincome_1990<14000
	replace hhfaminc_1990 = 13 if hhincome_1990>=14000 & hhincome_1990<15000
	replace hhfaminc_1990 = 14 if hhincome_1990>=15000 & hhincome_1990<17500
	replace hhfaminc_1990 = 15 if hhincome_1990>=17500 & hhincome_1990<20000
	replace hhfaminc_1990 = 16 if hhincome_1990>=20000 & hhincome_1990<22500
	replace hhfaminc_1990 = 17 if hhincome_1990>=22500 & hhincome_1990<25000
	replace hhfaminc_1990 = 18 if hhincome_1990>=25000 & hhincome_1990<27500
	replace hhfaminc_1990 = 19 if hhincome_1990>=27500 & hhincome_1990<30000
	replace hhfaminc_1990 = 20 if hhincome_1990>=30000 & hhincome_1990<32500
	replace hhfaminc_1990 = 21  if hhincome_1990>=32500 & hhincome_1990<35000
	replace hhfaminc_1990 = 22 if hhincome_1990>=35000 & hhincome_1990<40000
	replace hhfaminc_1990 = 23 if hhincome_1990>=40000 & hhincome_1990<50000
	replace hhfaminc_1990 = 24 if hhincome_1990>=50000 & hhincome_1990<75000
	replace hhfaminc_1990 = 25 if hhincome_1990>=75000 & hhincome_1990!=.

keep if year==`i'

collapse (sum) asecfwt, by(hhfaminc_1990 division)
rename asecfwt hh_count

egen sumhh=total(hh_count)
gen weight=hh_count/sumhh
save `i'_cps_division,replace
}



*** predictefd kwh

clear
forval i=1990/2020  {
clear
use recs_1990
capture drop if kwh==999999

bysort division moneypy: egen sumhh=total(nweight)
bysort division moneypy: gen weight=nweight/sumhh
bysort division moneypy: egen meankwh=total(weight*kwh)

 collapse meankwh,by(moneypy division)

quietly drop if division==. 
rename moneypy hhfaminc_1990
sort hhfaminc_1990 division
quietly destring hhfaminc_1990,replace

quietly merge 1:1 hhfaminc_1990 division using  `i'_cps_division

quietly drop if meankwh==.
quietly drop if hh_count==.

quietly sum meankwh[aweight=weight]
quietly return list
dis `i',r(mean), r(sd)

}



***==== all variables figure a5***===


clear
use income_b1990.dta

forvalues i=1990/2020{
use income_b1990.dta

* create real income groups
	gen hhfaminc_1990 = 1 if hhincome_1990<3000
	replace hhfaminc_1990 = 2 if hhincome_1990>=3000 & hhincome_1990<4000
	replace hhfaminc_1990 = 3 if hhincome_1990>=4000 & hhincome_1990<5000
	replace hhfaminc_1990 = 4 if hhincome_1990>=5000 & hhincome_1990<6000
	replace hhfaminc_1990 = 5 if hhincome_1990>=6000 & hhincome_1990<7500
	replace hhfaminc_1990 = 7 if hhincome_1990>=7500 & hhincome_1990<9000
	replace hhfaminc_1990 = 8 if hhincome_1990>=9000 & hhincome_1990<10000
	replace hhfaminc_1990 = 9 if hhincome_1990>=10000 & hhincome_1990<11000
	replace hhfaminc_1990 = 10 if hhincome_1990>=11000 & hhincome_1990<12500
	replace hhfaminc_1990 = 12 if hhincome_1990>=12500 & hhincome_1990<14000
	replace hhfaminc_1990 = 13 if hhincome_1990>=14000 & hhincome_1990<15000
	replace hhfaminc_1990 = 14 if hhincome_1990>=15000 & hhincome_1990<17500
	replace hhfaminc_1990 = 15 if hhincome_1990>=17500 & hhincome_1990<20000
	replace hhfaminc_1990 = 16 if hhincome_1990>=20000 & hhincome_1990<22500
	replace hhfaminc_1990 = 17 if hhincome_1990>=22500 & hhincome_1990<25000
	replace hhfaminc_1990 = 18 if hhincome_1990>=25000 & hhincome_1990<27500
	replace hhfaminc_1990 = 19 if hhincome_1990>=27500 & hhincome_1990<30000
	replace hhfaminc_1990 = 20 if hhincome_1990>=30000 & hhincome_1990<32500
	replace hhfaminc_1990 = 21  if hhincome_1990>=32500 & hhincome_1990<35000
	replace hhfaminc_1990 = 22 if hhincome_1990>=35000 & hhincome_1990<40000
	replace hhfaminc_1990 = 23 if hhincome_1990>=40000 & hhincome_1990<50000
	replace hhfaminc_1990 = 24 if hhincome_1990>=50000 & hhincome_1990<75000
	replace hhfaminc_1990 = 25 if hhincome_1990>=75000 & hhincome_1990!=.

keep if year==`i'

xtile ageg = age, nq(3)
xtile famsizeg=famsize,nq(3) 
gen white=(race==100)

collapse (sum) asecfwt, by(hhfaminc_1990 division famsizeg ageg white)
rename asecfwt hh_count

egen sumhh=total(hh_count)
gen weight=hh_count/sumhh
save `i'_all,replace
}

**check size
forval i=1990/2020 {
clear
use `i'_all.dta
di _N
}
*** predictefd kwh

clear
forval i=1990/2020  {
clear
use recs_1990

xtile ageg = hhage, nq(3)
xtile famsizeg=famsize,nq(3) 
gen white=(origin==1)

bysort division moneypy famsizeg ageg white: egen sumhh=total(nweight)
bysort division moneypy famsizeg ageg white: gen weight=nweight/sumhh
bysort division moneypy famsizeg ageg white: egen meankwh=total(weight*kwh)
collapse meankwh,by(moneypy division ageg famsizeg white)


quietly drop if division==. 
rename moneypy hhfaminc_1990
sort hhfaminc_1990 division
quietly destring hhfaminc_1990,replace

quietly merge 1:1 hhfaminc_1990 division famsizeg ageg white using  `i'_all
quietly drop if meankwh==.
quietly drop if hh_count==.

*sum kwh // this sum is not estimated/ estimated is by weights, this is just summary
quietly sum meankwh[aweight=weight]
quietly return list
dis r(mean)

}




